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We present a numerical analysis of the effect of particle elongation on the quasistatic behavior 
of sheared granular media by means of the Contact Dynamics method. The particle shapes are 
rounded-cap rectangles characterized by their elongation. The macroscopic and microstructural 
properties of several packings subjected to biaxial compression are analyzed as a function of particle 
elongation. We find that the shear strength is an increasing linear function of elongation. Performing 
an additive decomposition of the stress tensor based on a harmonic approximation of the angular 
dependence of branch vectors, contact normals and forces, we show that the increasing mobilization 
of friction force and the associated anisotropy are key effects of particle elongation. These effects are 
correlated with partial nematic ordering of the particles which tend to be oriented perpendicular to 
the major principal stress direction and form side-to-side contacts. However, the force transmission 
is found to be mainly guided by cap-to-side contacts, which represent the largest fraction of contacts 
for the most elongated particles. Another interesting finding is that, in contrast to shear strength, 
the solid fraction first increases with particle elongation, but declines as the particles become more 
elongated, ft is also remarkable that the coordination number does not follow this trend so that the 
packings of more elongated particles are looser but more strongly connected. 

PACS numbers: 45.70.-n,83.80.Fg,61.43.-j 



I. INTRODUCTION 



Since a few years, the research for a better understand- 
ing of the complex rheology of granular media is enriched 
by an increasing focus on nonspherical particles [U-Q- 
The wide-spread use of spherical or disk-like particles has 
been motivated by the fact that the rheology of granular 
media is basically governed by the collective contact in- 
teractions of the particles so that the particle shape can 
be viewed as a secondary effect. In practice, both in ex- 
periments and discrete element simulations, the spherical 
or circular particles, such as glass beads and disks, are 
easier to handle and the results are generally more di- 
rectly amenable to theoretical analysis. However, owing 
to the fast progress in experimental and numerical tech- 
niques during the last decade, there is now a wide scope 
for the investigation of materials composed of more com- 
plex particle shapes. In this respect, the model granular 
media with spherical particles provide a reference ma- 
terial for understanding the rheology when the particle 
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shapes deviate from a spherical or circular shape 

A wide variety of particle shapes can be found in na- 
ture and industry: elongated and platy shapes, e.g. in 
biomaterials and pharmaceutical applications, angular 
and facetted shapes, e.g. in geomaterials, and nonconvex 
shapes, e.g. in sintered powders. The behavior under 
various types of loading is strongly influenced by parti- 
cle shape. Rounded particles enhance flowability whereas 
angular shape is susceptible to improve shear strength. 
In many applications, the particle shapes need to be op- 
timized in order to increase performance [l2| - tl7j . These 
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trends are generally explained in qualitative terms and 
linked with the jamming of the particles. 

The effect of particle shape is mediated by the spe- 
cific granular texture (or fabric) induced by each particle 
shape. For example, it is found that hard ellipses can be 
jammed even though they are underconstrained (l8l - [23j . 
In general, the anisometric or elongated particle shapes, 
such as spheroids and sphero-cylinders, tend to develop 
orientational order affecting force transmission and fric- 
tional behavior 0, 124H26I ] . This "nematic" ordering oc- 
curs while, in contrast to liquid crystals, the particles 
interact only via contact and friction [27| . 

In a sheared granular material, the local equilibrium 
structures are generically anisotropic in terms of contact 
directions and forces I28l - l34j . It was recently shown 
that the fabric anisotropy in a sheared granular assembly 
crucially depends on particle shape [HI [H| • In the case 
of polygonal and polyhedral particles, due to large con- 
tact area of side-to-side contacts, the fabric anisotropy 
appears to be marginal compared to force anisotropy 
[pU [lT| . Those contacts play a major role in force trans- 
mission by accommodating long force chains that are ba- 
sically unstable in a packing composed of spheres. 

The force and fabric anisotropics are at the origin of 
the enhanced shear strength of materials composed of 
nonspherical particles [ll|, HH HH HH • The particle shape 
affects the compactness and dilatancy of granular mate- 
rials. A nontrivial effect, evidenced recently by exper- 
iments and numerical simulations for spheroids, is the 
finding that the solid fraction is not a monotonous func- 
tion of the aspect ratio [l9l - [22l |36| . The solid fraction 
increases linearly to a maximum and then declines in 
inverse proportion to the aspect ratio [37j • In powder 
processing, the particle shape appears also to be an im- 
portant parameter controlling the flowability, discharge 
mDatpa.and compaction of powders [38l [39j| . 
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In this paper, we use contact dynamics simulations 
to investigate the rheology of large packings of elon- 
gated particles with increasing aspect ratio. The particles 
are rectangles with rounded caps to which we will refer 
as Rounded-Cap Rectanglular (RCR) particles. These 
particles may be considered as 2D analog of sphero- 
cynlinders. The RCR shape can be characterized by a 
single aspect ratio a or, as we shall see, by an elongation 
parameter r\ varying from to 1 as the particle shape 
varies continuously from a circle to an thin line. We are 
interested both in the properties of the static packings of 
RCR particles prepared by isotropic compaction without 
friction and in the stress-strain behavior under biaxial 
compression with finite friction between particles. 

The macroscopic behavior is studied in terms of the 
internal angle of friction and solid fraction for different 
values of r\. We find a nonmonotonous variation of the 
solid fraction and a nearly linear increase of the inter- 
nal angle of friction with ry. In order to understand the 
origins of this behavior and the role of particle shape, 
we perform a detailed analysis of the microstructure and 
stress transmission. We consider the organization of the 
particles and contacts in the simulated packings, as well 
as the stress transmission by means of a harmonic repre- 
sentation of the stress tensor in terms of force and fabric 
anisotropics. The microstructure is increasingly domi- 
nated by a short-range nematic ordering of particle ori- 
entations as i) increases. We show that the internal angle 
of friction is influenced by this ordering via an increasing 
anisotropy of friction forces and contact orientations with 
the elongation parameter. For all values of the latter, the 
harmonic approximation provides an excellent fit to the 
shear stress. 

An important feature of RCR particles is that, like 
polygonal particles, they have lineal edges and can thus 
form side-to-side contacts as well as side-to-cap and cap- 
to-cap contacts. Hence, in a packing of RCR particles, 
the texture can characterized by the networks of these 
various contact types, and the influence of the shape pa- 
rameter on force transmission and shear strength may be 
analyzed in terms of these contacts and more specially 
the side-to-side contacts which are expected to play a 
stabilizing role in the packing. 

In the following, we first introduce our numerical ap- 
proach in Section [TTJ Then, in Section IIII1 the stress- 
strain behavior is presented for different values of r\. The 
microstructure is analyzed in Section HVl in terms of con- 
nectivity, orientations of the particles and the contact 
network. We also introduce the harmonic approximation 
of the stress tensor allowing us to track the origins of the 
internal angle of friction via force and fabric anisotropics. 
In Section |V1 we present an additive decomposition of 
the connectivity, anisotropics and forces as a function 
of different contact types. The force distributions are 
presented in Section [VIl In section Ivm we analyze the 
structure of force networks with cap-to-cap, cap-to-side 
and side-to-side contacts. We conclude with a summary 
and discussion of the most salient results of this work. 



II. SYSTEM DESCRIPTION AND NUMERICAL 
PROCEDURES 

In this section, we briefly introduce the contact dy- 
namics (CD) method used for the simulations and the 
numerical procedures for sample preparation. 

A. Contact Dynamic method 

The simulations were carried out by means of the con- 
tact dynamics (CD) method @, 113411 ■ The CD method 
is a discrete element approach based on a nonsmooth ap- 
proach in which an integrated form of the equations of dy- 
namics. The integration time interval corresponds to the 
time step and may involve discontinuous variation of the 
velocities due to collisions. The frictional and collisional 
interactions are described as complementarity relations 
between the relative velocities between particles and the 
corresponding momenta at the contact points without 
elastic or viscous regularization. Thus, the condition of 
geometrical contact between two particles is expressed by 
the following mutually exclusive alternatives: 

/„ ^ and u n = 0, , , 

f n = and m„>0. 1 ' 

where f n is the normal contact force and u n the rela- 
tive normal velocity between two particles in contact is 
counted positive when they move away from each other. 

In the same way, the Coulomb friction law involves 
three mutually exclusive conditions: 

ft = ~nf n and u t > 0, 
~M/n ^ ft < M/n and u t = 0, (2) 
ft = M/n and u t < 0, 

where u t is the sliding velocity at the contact, fi is the 
friction coefficient and ft is the friction force. Remark 
that this relation cannot be reduced to a (mono)valucd 
functional dependence between the two variables as as- 
sumed in the Molecular Dynamics (MD) method. 

The above formulation is implicit in the sense that the 
complementarity relations should be satisfied for the ve- 
locities at the end of each time step. An iterative algo- 
rithm based on a nonlinear Gauss-Seidel scheme is used 
to solve the system of equations and complementarity 
relations for contact forces and particles velocities. The 
uniqueness is not a priori guaranteed for perfectly rigid 
particles. However, by initializing each step of calcula- 
tion with the forces calculated in the preceding step, the 
set of admissible solutions shrinks to fluctuations which 
are basically below the numerical resolution. 

The CD method is particularly suitable for the sim- 
ulation of rigid nearly undeformable particles. In this 
limit, the MD method requires steep interaction poten- 
tials and thus very small time steps. Nevertheless, several 
comparisons between the two methods suggest that both 
methods are equally valid and efficient for the simulation 
of granular materials @, H3, EH . 



FIG. 1: Shape of a Rounded-Cap Rectangle (RCR). 

B. Simulation of RCR particles 

We model the RCR particle as a juxtaposition of two 
half-disks of radius R' with one rectangle of length L 
and width 2R'; see Fig. [TJ The shape of a RCR parti- 
cle is a circle of radius R' for L = 0. The aspect ratio 
a — (L + 2R') / (2R 1 ) is 1 in this limit and increases with 
L for a fixed value of R' . In this paper, we use an alter- 
native parameter describing the deviation of the particle 
shape from a circle. Let R be the radius of the circle 
circumscribing the particle. We have R = L/2 + R'. The 
radius R is also that of the inscribed circle. Hence, the 
deviation from a circular shape can be characterized by 
AR, = R — R' = L/2. We use the dimensionless parame- 
ter -q defined by 



It varies from 77 = 0, for a circle, to 1 corresponding to a 
line. We will refer to r\ as the elongation parameter as in 
rock mechanics (Hoj . 

For the detection of the contacts between two RCR 
particles, we use the schema shown in Fig. [2j Three types 
of contact can be distinguished: cap-to-cap (cc), cap-to- 
side (cs) and side-to-side (ss). The contacts between 
the particles are thus detected separately for the pairs of 
circles and rectangles. In general, in the CD method ss 
contact between two rectangles is treated as composed 
of two point contacts and the contact laws (HJ and © 
are applied separately to each point. The choice of these 
points does not affect the resultant force and its point 
of application. Hence, for RCR particles, as shown in 
Fig. [2] ss contact is composed of four point contacts : 
two points due to the rectangle-rectange interface and 
two points due to the cc contacts. Thus, four forces are 
calculated by the CD algorithm but only their resultant 
and application point are material. 

The detection of line contacts between rectangles 
was implemented through the so-called shadow overlap 
method devised initially by Moreau [lol |4o| for polygons. 
The reliability and robustness of this method have been 
tested in several yea rs of previous applications to gran- 
ular materials [3I llOL [Til O |H, [5ll [52j]. This detection 
procedure is fairly rapid and allows us to simulate large 
samples composed of RCR particles. For our simulations, 
we used the LMGC90 which is a multipurpose software 
developed in Montpellier, capable of modeling a collec- 
tion of deformable or undeformable particles of various 



FIG. 2: Representation of cap-to-cap, cap-to-side and side- 
to-side contact and they will be referred as cc contacts, cs 
contacts and ss contact, respectiveley. 



shapes (spherical, polyhedral, or polygonal) by means of 
the contact dynamics (CD) method |45| . 



C. Sample preparation 

We prepared 8 different packings of 13000 RCR par- 
ticles with r) varying from to 0.7 by steps of 0.1. The 
radius R of the circumscribing circle defines the size of 
a RCR particle. In order to avoid long-range ordering in 
the limit of small values of 77, we introduce a size poly- 
dispersity by taking R in the range [R m in, Rmax] with 
Rmax = %Rmin with a uniform distribution in particle 
volume fractions. 

All samples are prepared according to the same pro- 
tocol. A dense packing composed of disks (77 = 0) is 
first constructed by means of a layer- by- layer deposition 
model based on simple geometrical rules [5314561 ] . The 
particles are deposited sequentially on a substrate. Each 
new particle is placed at the lowest possible position at 
the free surface as a function of its diameter. This proce- 
dure leads to a random close packing in which each parti- 
cle is supported by two underlying particles and supports 
one or two other particles. For rj > 0, the same packing 
is used with each disk serving as the circumscribing circle 
of a RCR particle. The RCR particle is inscribed with 
the given value of 77 and random orientation in the disk. 

Following this geometrical process, the packing is com- 
pacted by isotropic compression inside a rectangular 
frame of dimensions Iq x ho in which the left and bottom 
walls are fixed, and the right and top walls are subjected 
to a compressive stress Co- The gravity g and friction 
coefficients fi between particles and with the walls are 
set to zero during the compression in order to avoid force 
gradients and obtain isotropic dense packings. Fig. [3] 
displays snapshots of the packings for several values of 77 
at the end of isotropic compaction. 

The isotropic samples are then subjected to vertical 
compression by downward displacement of the top wall 
at a constant velocity v y for a constant confining stress 
<To acting on the lateral walls. The friction coefficient \x 
between particles is set to 0.5 and to zero with the walls. 
The simulations were run with a time step of 2.10~ 4 s. 
The CPU time was 5.10~ 4 s per particle and per time 
step on an AMD processor. Since we are interested in 
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FIG. 3: Examples of the generated packings at the initial 
state. 

quasistatic behavior, the shear rate should be such that 
the kinetic energy supplied by shearing is negligible com- 
pared to the static pressure. This can be formulated in 
terms of an inertia parameter I defined by [57j : 

(fl 



where e = y/y is the strain rate, m is the particle mass 
and p is the mean pressure. The quasistatic limit is char- 
acterized by the condition In our simulations, I 
was below 10~ 3 . 



III. STRENGTH AND DILATANCY 

In this section, we consider the stress-strain and 
volume-change behavior as a function of the shape pa- 
rameter ij. We need to evaluate the stress tensor and solid 
fraction during deformation from the simulation data. 
For the stress tensor, we start with the tensorial moment 
M % of each particle i that is defined by [H, m : 

where /° is the a component of the force exerted on 
particle i at the contact c, is the /3 component of the 
position vector of the same contact c, and the summation 
runs over all contact neighbors of particle i (noted briefly 
by c £ i). The average stress tensor <x in the volume 
V of the granular assembly is given by the sum of the 
tensorial moments of individual particles divided by the 
volume [HHIi : 

where t c is the branch vector joining the centers of the 
two touching particles at the contact point c. Remark 
that the first summation runs over all particles whereas 
the second summation involves the contacts, each contact 
appearing only once. 

Under biaxial conditions with vertical compression, we 
have <j\ > cr 2 , where the a a are the stress principal val- 
ues. The mean stress p and stress deviator q are defined 
by: 

P = ^(0-1+02), (7) 

Q = - °2)- (8) 

For our system of perfectly rigid particles, the stress state 
is characterized by the mean stress p and normalized 
shear stress q/p. 

The strain parameters are the cumulative vertical, hor- 
izontal and shear strains e 1; e 2 and e q , respectively. By 
definition, we have 



[ h dh' , / 













where ho is the initial height and Ah = ho — h is the total 
downward displacement, and 

f l dV f Al\ 

Ej = i„- = '"( 1+ v)' (10) 
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where Iq is the initial box width and Al — I — Iq is the 
total change of the box width. The cumulative shear 
strain is then defined by 

e q = e 1 -e 2 - (11) 
Finally, the cumulative volumetric strain e p is given by 



£l + £2 



f v dV . ( 













(121 



where Vq — l^h^ is the initial volume and Av = v — vq is 
the cumulative change of solid fraction. 

Figure 2] shows the normalized shear stress q/p as a 
function of shear strain e q for different values of 77. The 
jump observed at e q = reflects both the rigidity of the 
particles and high initial solid fraction of the samples 
(see below). In all cases, the shear stress passes by a 
peak before relaxing to a stress plateau corresponding to 
the so-called "residual state" in soil mechanics (60[. We 
remark that the residual shear stress increases with 77. 

The internal angle of friction <p* , representing the shear 
strength of the material, is defined from the mean value 
(q/ p)* of the normalized shear stress in the residual state 
by [60] 



sin ip 



(I)" 



(13) 



Fig. @] shows the variation of sin tp* as a function of a 
and 77. We see that the shear strength is an increasing 
nonlinear function of the aspect ratio, but, interestingly, 
it varies linearly when plotted versus the elongation pa- 
rameter. Hence, we have 



sin p* = sin p^ + k ij = sin Pq + k 



1 

1 

a 



(14) 



This observation indicates that the evolution of shear 
strength reflects more directly shape elongation than as- 
pect ratio. In the following, we will use 77 as shape pa- 
rameter. 

Figure [5] (a) displays the cumulative volumetric strain 
e p as a function of e q for different values of 77. Start- 
ing with an initially dense state, all packings dilate and 
hence the volume increases. For r\ < 0.4, a plateau is 
reached beyond e q — 0.3, corresponding to a state of 
isochoric deformation. For larger aspect ratios, the di- 
latation continues even at very large deformations. This 
is an indication of an inhomogeneous dilation due to the 
formation of shear bands in the bulk, which is enhanced 
by the elongated shape of the particles. Since different 
parts of the packing undergo differential volume change, 
longer shearing is required to reach a fully dilated state 
for the whole packing. The initiation and evolution of 
shear bands for different particle elongations will be re- 
ported in more detail elsewhere. 

Figure [5] (b) displays the solid fraction v as a func- 
tion of 77 at different levels of shear deformation e q . It 
is remarkable that, at all levels of deformation, the solid 
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FIG. 4: Normalized shear stress q/p as a function of cumula- 
tive shear strain e q for different values of the shape parameter 
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FIG. 5: Internal angle of friction p* as a function of aspect 
ratio a (a) and elongation rj. The error bars represent the 
standard deviation in the residual state. 



fraction increases with 77, reaches a maximum at r\ ~ 0.4 
and then declines as 77 further increases. We note that 
solid fractions as large as 0.90 are reached for 77 = 0.4 in 
the initial state. A similar nonmonotonous behavior was 
observed for packings of ellipses or ellipsoidal particles 
[l9l [20I [22j . This is somewhat a counterintuitive finding 
as the shear strength (a monotonous function of 77) does 
not follow the trend of solid fraction (nonmonotonous). 
This behavior is clearly not related to shear localization 
since it is observed at all levels of deformation includ- 
ing the initial isotropic state v$ = v(e q — 0) where the 
packings are homogeneous. 

A rapid fall-off of solid fraction for elongated parti- 
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FIG. 7: Dilatancy angle ip versus internal angle of friction p> 
for different values of n. 



In the following, we focus on the microstructure and 
force transmission that provide a key to a better under- 
standing of the physical mechanisms underlying the effect 
of particle shape on the shear strength. 



IV. GRANULAR TEXTURES 



FIG. 6: Cumulative volumetric strain e v as a function of shear 
strain s q (a); Solid fraction as a function of particle shape 
parameter n (b) at different levels of shear strain. 



cles in 3D was observed at very large aspect ratios and 
attributed to the excluded volume due to disorder, pre- 
dicting a fall-off in inverse proportion to the aspect ratio 
[37l Irlll 62] . The initial rapid increase of solid fraction, as 
observed in Fig. [51 reveals that excluded-volume effects 
are not the prevailing mechanism at low aspect ratios. 
In this limit, slight deviations from spherical shape have 
strong space-filling effect on the packing although the ex- 
cluded volume increases at the same time and becomes 
dominant at very large aspect ratios. 

The volumetric deformation can also be expressed in 
terms of the dilatancy angle ip defined by 63] : 



sin-0 



de p 
de„ 



(15) 



The plot of ip as a function of ip, the so-called "stress- 
dilatancy diagram" , is shown in Fig. [7]for different values 
of rj. We observe a linear correlation between ip and ip 
irrespective of the value of r\. We have 



ip ~ ip* + ip. 



(16) 



This is a particularly simple relation compared to several 
models proposed in soil mechanics [HI, l64j . It reflects 
the "non-associated" character of granular plasticity, an 
associated behavior implying simply gp = ip, which is un- 
realistic for granular materials [63l . 16514671 ] . According to 
relation (fT6)l . the dilatancy angle vanishes in the residual 
state where (p = ip*. Recent work on cohesive and granu- 
lar packings of polygonal particles in 2D is in agreement 
with this correlation \3U, mm ■ 



In this section, we investigate the general organization 
(texture) of our packings of RCR particles in terms of 
particle orientations and contact network. This will allow 
us to quantify the effect of the elongation parameter and 
its connection with shear strength. 



A. Particle orientations 

The principal feature of elongated particles is their ori- 
entational degree of freedom. The particle orientation is 
represented by a unit vector m as shown in the inset 
to Fig. [9] In 2D, it is parametrized by a single an- 
gle i?. Let T> ($) be the set of particles with direction 
i? 6 [$ — S'd/2; d + 8d/2] for angle increments Sd, and 
Np($) its cardinal. The probability P^(^) of the orienta- 
tions of particles is given by 



(17) 



where N p is the total number of particles. 

Figure [S] displays a polar representation of P§ (■&) for 
rj = 0.7 at various stages of deformation. In the initial 
state, corresponding to an isotropic stress state, the dis- 
tribution is anisotropic with privileged direction i} p close 
to 7r/2. This means that for elongated particles, the par- 
ticle orientations are not fully correlated with the stress 
state so that the resulting particle orientation anisotropy 
depends on details of the assembling process that can 
not be controlled by simply subjecting the particles to 
isotropic stresses from the boundary. 

The priviliged direction rotates as a result of vertical 
compression and becomes horizontal (parallel to the mi- 
nor principal stress direction) in the residual state. The 
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FIG. 9: Particle orientation anisotropy dp JUS £L function of 
shape parameter tj at different stages of shear e q . 



FIG. 8: Polar representation of the probability density P$ 
of particle orientations i9 for rj = 0.7 at various stages of 
deformation e q . The symbols are the simulations data. Solid 
lines represent harmonic fit according to equation (|18[) . 



distribution are nicely fitted by harmonic approximation 
corresponding to the lowest order terms of the Fourier 
expansion of P#($) [H, Iz3| : 



W) - 7^(1 + ^008(2(1?-^)), 
Z7T 



(18) 



where a p represents the anisotropy of the distribution 
and 6 a is the major principal stress direction. The choice 
of 6 a as reference direction is motivated by the observa- 
tion that the privileged orientation of the particles tends 
to align itself with the minor principal stress direction. 
Hence, negative values of a p mean that the particles are 
preferentially oriented perpendicular to the major prin- 
cipal stress. 

We plot the particle orientation anisotropy a p in Fig. 
IHlas a function of rj at different stages of shear e q . We see 
that the particle orientations are isotropic for 77 < 0.4 at 
the initial state and they become increasingly anisotropic 
as r\ increases beyond 0.4. At nearly all stages of shear, 
a p is negative, and at most advanced stage, i.e. corre- 
sponding approximately to the residual state, it becomes 
nearly independent of 77. The large absolute value of a p 
(~ 0.35) suggests that many particles are aligned in hor- 
izontal layers just as in nematic order. One example of 
this nematic ordering is shown in Fig. [TU] for two differ- 
ent values of 77. This ordering may be attributed to the 
favored mechanical equilibrium of the particles under the 
action of vertical stress and enhanced by boundary align- 
ment of the elongated particles 0, HI, 0, [7l[ . This point 
will be analyzed more deeply below. 



B. Particle connectivity 

The primary statistical quantity describing the contact 
network is the coordination number z (average number 



of contacts per particle). For our elongated particles, 
each side-to-side contact is counted as one contact even 
if side-to-side contacts are treated as four point contacts 
belonging to the contact segment (see section [TTJ) . The 
floating particles with no force-bearing contact (i.e. with 
less than two active contacts) are thus removed from the 
statistics. The fraction of floating particles decreases lin- 
early with 77 from 17% for 77 = to 10% for 77 = 0.7. 

Figure [TT1 displays the evolution of z as a function of 77, 
in the initial and residual states. The initial-state value 
z = 4 corresponds to a frictionless packing of circular 
particles in the isostatic state with z = 2d where d is 
space dimension (indeed, the packings where prepared 
by setting the friction coefficient to zero). However, as 77 
increases, z increases to a plateau value of ~ 5.6. This 
is in agreement with recent work showing that large dis- 
ordered jammed packings are isostatic only for disks or 
spheres l -> 2,'J. \72\ . For nonspherical or noncircular par- 
ticles, we have z < d(d + 1). 

Numerical results for frictional or frictionless systems 
of rigid or deformable disks and sphere s Il8l . [7^h77|. as 
well as with ellipses and spheroids [HW23I ] confirm this 
point. In the residual state, the mean value of z is below 
that in the isotropic state, and it grows from 3 to 5 with 
77. It is interesting to note that z does not follow the solid 
fraction which, as we have seen before, is nonmonotonous 
as a function of 77. This means that for large aspect ratios, 
the packings are loose but well connected. 

The connectivity of the contact network can be char- 
acterized in more detail by the fraction P(c) of particles 
with exactly c contact neighbors. Fig. [T2l shows P{c) in 
the residual state for different values of 77. The distri- 
bution is increasingly broader as 77 becomes larger. The 
particles can have as many as 10 contact neighbors at 
77 = 0.7. This is allowed both by the geometry and poly- 
dispersity of the particles as shown by a typical grey-level 
map of particle connectivities in Fig. 1131 For 77 = 0, we 
observe a peak at c = 3. This peak slides gradually to 
c = 5at77 = 0.7as observed also for z, which is, by 
definition, the mean value of c for force-bearing contacts: 
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FIG. 10: (Color online) Color level map of particle orienta- 
tions for n — 0.2 (up) and n = 0.7 (down) in the residual 
state. 



C. Force and texture anisotropies 

Equation ^ shows that the expression of stress tensor 
is an arithmetic mean involving the branch vectors £ and 
contact force vectors /. This means that for the analysis 
of stress transmission and shear strength from a particle- 
scale viewpoint we need a statistical description of these 
quantities. 

A common approach used by various authors is to ex- 
press branch vectors and contact force orientations in 
terms of the contact direction, i.e. in the local con- 
tact frame (n,t), where n is the unit vector perpendic- 
ular to the contact plane, and t is an orthonormal unit 
vector oriented along the tangential force [3, H, H3, [M- 
M, HI, [M [ll-lil ; see figure [HJa). The components of 
the branch vector and contact force are expressed in this 
frame as: 



£ n n + £ t t, 
f n n + f t t, 



(19) 



N 



» Initial state 
■ Critical state 



I I I I I I I I I I I 

-Ul 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 



FIG. 11: Initial and residual coordination numbers as a func- 
tion of shape parameter n. The error bars represent the stan- 
dard deviation in the residual state. 




FIG. 12: Connectivity diagram for each samples expressing 
the fraction P(c) of particles with exactly c contacts in the 
residual state. Note that the floating particles (i.e. with more 
than one active) are removed from the statistics 



nents of the branch vectors, and f n and ft the normal 
and tangential components of the contact force. 

Remark that only for disks or spherical particles we 
have t = in where £ is the length of the branch vec- 
tor. A consequence of noncircular or nonspherical par- 
ticle shape is to dissociate the contact frame from the 
branch vector frame (n',i'), where n' is the unit vector 
along the branch i and t' is the orthoradial unit vector 
[UHIl ; see figure IT3rb). We express the components of 
the branch vector and contact force also in this frame: 



Z n >n, 

f n ,n' + f t ,t', 



(20) 



where £ n and It are the normal and tangential compo- 



where f n > and /e are the radial and orthoradial compo- 
nents of the contact force, and l n i = I. 

In two dimensions, let 9 and 9' be the orientations 
of of n and n' , respectively. From the numerical data, 
we can evaluate the probability density functions Pe{9) 
and Pgr (6') of contact and branch vector orientations, 
respectively, as well as the angular averages of the force 
components (/„)(#), (/t)(0), (/n'X^). and {fv){9') and 
branch vector components (£ n }(9), (£ t )(9), (£„/)(#'). In 
the absence of an intrinsic polarity for n and n', all these 
functions are 7r-periodic. The insets to Figs. [TSl HU and 
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FIG. 13: (Color online) Color map of particle connectivities. 
Color intensity is proportional to coordination number. 





FIG. 14: Contact frame 
(n',f) (b) 



(n, t) (a) and intercenter frame 



[T7]display polar representations of these functions for 77 = 
0.7 at the end of shearing. All these angular functions 
are generically anisotropic. The peak values occur along 
the axis of compression (8 = tt/2) for Pg, Pg> , (/„) and 
(/„/}, and along the axis of extension (8 = 0) for (l n ) 
and (£„,')■ The maxima for the tangential components 
occur in the direction of 7r/4 with respect to the axis of 
compression. 

The simple shapes of the above functions suggest that 
they can be approximated by their Fourier expansions up 
to the second term Fiji [H, IS " 



f Pe(6) 

<M(0) 
<M(e) 
</n»)(e) 
I </t->(e) 



£={1 + <? C cos 2(6 -0*)}, 
{£){! + ain* cos2(e-e in .)}, 
(£)a lt *sm2(Q-e lt ,), (21) 
(/){l + o /n . cos2(9-e /n *)} 
{f)a ft . sin2(e-e /t .), 



where stands either for 8 or for 8' depending on the lo- 
cal frame used. The (£) is mean length of branch vectors 
and (/) is the mean force, (a*, a;„* , aw , , aft*) — 
(a c ,ai n ,ait,af n ,af t ) and (6*, ©;„», ©it*> ©,/V, ©ft*) = 
{0c,0i n ,6it,6f n ,8 ft) are the anisotropy parameters and 
the angle of privileged direction of each function 
in the frame (n,i). In the same way, we have 



(a' c ,ain',aitt,af n t,a,ft>) and 
= (9' c ,8in' ,9if ,9f n ' ,8ff) in 



(a*> a in*> ajt*» o/t»*j o/t*) : 

(e;,e In .,e,t.,e/„.,e/ t . 

the frame (n', t'). 

Note that, by construction, we have aw — and 
= 0. In the following, we will refer to a c as 
contact anisotropy, to a' c as branch vector orientation 
anisotropy, to (ai n * , aw ) as branch length anisotropics 
and to (a /„* , a/t* ) as normal and tangential or radial 
and orthoradial force anisotropics depending on the lo- 
cal frame 0, |35[. These harmonic approximations are 
well fit to our data as shown in the insets to Figs. [T5l [T6l 
and El 

In practice, it is convenient to estimate the above 
anisotropics from the following fabric and force tensors 



F, 



aft 



in 
A a /3 

It* 
Xa/3 

fn* 
X a /3 

X a /3 



±Jn* a n* P e (Q)dO, 


ijV„*)(©)<n£Pe(©)de, 


^](£t*K®)n* a t* P e (O)dO, 


^j(fn*)(&)n* a n* p P (e)de, 



(22) 



where a and (3 design the components in the considered 
frame. Note that, by construction, we have XaB = 
From equations ([21]) and (|22|) , assuming that in a sheared 
state 8* = /n . = Off = CT , Q in * = K . = or 9 a , 
the following relations are easily obtained: 



ain* 
aw 
a fn - 
a ft * 



2(F 1 * - 
2(x!"* 

2(xi* 



2(X1 



Fi)/{F{ 

■x%)/(x l i 



ln ')/(x[ n ' 



Fi), 



X 2 
J" 



xt, 

/n *)/(x{"* J 



x l f 



Kx{ -xD/ixl + x 2 r ) 



xi n '" 



where % 



l* 



In* 



At* 



a* - 0/ n * , 
. ( 23 ) 

X"" + X" , X J = X J " + X ft and the 
indices 1 and 2 refer to the principal values of each tensor. 
By construction, we have (F* + F 2 * ) = 1, (xi +X2) = {fy 
and (x{ + X2 ) = (/)■ Note that a*, af n * and aft* are 
always positive whereas a;„« and aw are negative. We 
have 0;„. = and ®w — 0,. 

Figure [15] displays the variation of contact anisotropy 
a c and branch vector orientation anisotropy a' c , both av- 
eraged in the residual state, as a function of 77. We ob- 
serve two distinct behaviors: a c increases quickly from 
0.3 to 0.7 with 77 whereas, after a slight increase, a' c de- 
clines to nearly for 77 = 0.7. It is often admitted that 
the shear strength in granular media is a consequence 
of the buildup of an anisotropic geometrical structure 
due to friction and steric exclusions between particles 
[H S Izlj H3413|- But, here we have two different struc- 
tural anisotropics a c and a' c that vary in opposite direc- 
tions as 77 is increases. Hence, when the granular struc- 
ture is complex as in our packings of nonspherical particle 
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FIG. 15: Contact anisotropy a c and branch vector anisotropy 
(X c cLS cL function of shape parameter rj averaged in residual 
state. The error bars represent the standard deviation in the 
residual state. The inset shows the angular probability den- 
sities Pe(B) in black and P e i(8') in red for r\ — 0.7 calculated 
from the simulation data (points) together with the harmonic 
approximation (lines). 
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FIG. 16: Normal and tangential branch length anisotropics 
din and an and branch length anisotropy a in / as a function of 
shape parameter rj in the residual state. The error bars rep- 
resent the standard deviation in the residual state. The in- 
set shows the angular average functions (£i n )(6), (tu){0) and 
{lin>){@) m black, red and green, respectively, for rj = 0.7 cal- 
culated from the simulation data (points) and approximated 
by harmonic fits (lines). 



shapes, the choice of the statistical representation of the 
granular structure has to be specified [ill, HH ■ This point 
will be addressed in more detail in section I VIII 

The branch vector length anisotropies et/„, an and ai n >, 
averaged in the residual state, are plotted in FigOj)] as a 
function of rj. These parameters are negligibly small at 
small values of 77, i.e. for nearly circular particles, and 
decline to negative values as rj is increased. This means 
that the particles tend to form longer branch vectors with 
their neighbors in the direction of extension, suggesting 
that the particles touch preferentially along their minor 
axes when the contact orientation is close to the compres- 
sion axis, and along their major axis when the contact 
orientation is close to the extension axis; see section IVlIl 
It is also remarkable that a/„ ~ ai n / whereas |a;„| < |a/ t | 
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FIG. 17: Normal and tangential force anisotropies a„ and 
at and radial and orthoradial force anisotropies a' n and a' t as 
a function of rj in the residual state. The error bars repre- 
sent the standard deviation in the residual state. The inset 
shows the angular average functions (/n)(#) and (f n '}(d) in 
black and green, respectively, (a) and (ft)(9) and (/t'}(#) in 
red and blue, respectively, (b) for rj — 0.7 calculated from 
the simulation data (points) together with the harmonic ap- 
proximation (lines). The error bars represent the standard 
deviation in the residual state. 

particularly for 77 ^ 0.4. 

The normal and tangential force anisotropics aj n and 
aft are plotted in Fig. [17] as a function of 77, together 
with the radial and orthoradial force anisotropies af n i 
and aft', averaged in the residual state. In contrast to 
contact anisotropy, we see that a/„ and af n > grow to- 
gether slowly until 77 = 0.4, then a/„ remains nearly con- 
stant whereas af n > increases. On the other hand, the 
anisotropy aff of orthoradial forces grows much faster 
with rj than the anisotropy aj t of tangential forces. Re- 
markably, from rj > 0.4 the orthoradial force anisotropy is 
higher than radial force anisotropy (a/e > a/„') whereas, 
even if the tangential force anisotropy increases with 77, it 
is still below the normal force anisotropy (a/ t < aj n ) and 
remains always below the contact anisotropy (a/„ < a c ). 
In other words, the force anisotropy described in terms of 
branch vectors reflects more sensitively the effect of par- 
ticle shape elongation than in terms of contact normal 
vectors. We will see below that this behavior is related 
to the mobilization of friction (section [Vj) and contact 
types (section IVlip . 

V. GEOMETRICAL AND MECHANICAL 
ORIGINS OF SHEAR STRENGTH 

The stress tensor as formulated in Eq.© is a func- 
tion of discrete microscopic parameters attached to the 
contact network. For sufficiently large systems, the de- 
pendence of volume averages on individual discrete pa- 
rameters vanishes [H, [13, EE] and the discrete sums can 
be replaced by integrals. According to Eq. ©, we have 

V<T a p=Y,fZ*P= N c(f*t(>), ( 24 ) 



where N c is the total number of contacts. By writing the 
average on the right hand side in integral form, we get 



: [ faifi Ptfdf d£, 

Jo. 



(25) 



where Pgf is the joint probability density of forces and 
branch vectors, n c is the number density of contacts and 
fl is the integration domain in the space (£, f). 

The integral appearing in equation (f25j) can be reduced 
by integrating first with respect to the forces and branch 
vector lengths. Considering the components of the forces 
and branch vectors in one of the two local frames (n, t) 
or (n', t') and neglecting the branch vector- force correla- 
tions, we get d,[lH,|33,|33: 



^ = n c J {(e n *)(e) n* a (e) + (i t *)(e) t}(e)} 

o 

{</«*>(©) <(&) + (/**>(©) ^( )} p ( ) dQ - 

(26) 



The expression of the stress tensor by equation (|26| 
makes appear explicitly the average angular functions 
representing the fabric and force states. Using the har- 
monic approximation (I21[) , equation (|26p can be inte- 
grated with respect to space direction 6 and the stress 
invariants p and q extracted. Assuming that the stress 
tensor is coaxial with the fabric and force tensors (f22j). 
we get the following simple relations: 



P 



\{a c + a ln + a u + cif n + aft) in (n,t) 
\{a! c + aw + o/n' + aft') in (n' , t'). 



(27) 



The assumption of coaxiality is natural since, even if the 
preferential orientations of the forces and branch vectors 
are not fully correlated, we observe that shearing tends to 
align the contacts and forces with the principal directions 
of the stress tensor. 

Figure [TH] displays the residual-state value of the nor- 
malized shear stress (q/p)* as a function of rj calculated 
both directly from the simulation data and from equation 
(|27| separately for the two local frames by using the val- 
ues of various anisotropics estimated from the simulation 
data. As we see, for both local frames, equation (|2"7|) pro- 
vides an excellent fit to the data for all values of ij. Note, 
however, that the second expression in equation (|27| is 
more simple than the first expression (4 anisotropic pa- 
rameters vs 5 anisotropic parameters) and the resulting 
fit appears to be more accurate. 

The two equations (|27|) are interesting as they re- 
veal distinct origins of shear strength in terms of force 
and texture anisotropics with two different decomposi- 
tions. Various anisotropies do not contribute equally to 
the shear strength. The dominant anisotropies are tex- 
ture anisotropies for projection on contact frame since 
»e + ain + o-it > o! c + ai n > and force anisotropies for pro- 
jection on the branch vector frame since a/„' + aft' > 
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FIG. 18: Normalized shear stress (q/p)* in the residual state 
as a function of r\ together with two analytical expressions 
given by equation (|27|) . The error bars represent the standard 
deviation in the residual state. 



af n + aft- The fact that the texture anisotropy prevails 
in the contact frame may be attributed to its strong cor- 
relation with particle orientations. Geometrically, for a 
particle oriented along a direction more contacts can 
be formed with the flat side of the particle with normals 
oriented along i? + 7r/2 than with its rounded caps. This 
is consistent with the observation that the particle ori- 
entations are strongly anisotropic with an anisotropy a p 
of negative sign (preferred direction along the extension 
axis) and the contact normal anisotropy a c is positive 
(along the compression axis) and increases with aspect 
ratio; see section HVl 

In Fig. [TO] two maps of radial forces are shown for 
packings with 77 = 0.2 and 77 = 0.7, respectively. In the 
presence of long parallel sides, the strong force chains are 
more tortuous. Hence, the stability of such structures 
requires strong activation of tangential forces. Indeed, 
we remark that the orthoradial force anisotropy is above 
the radial force anisotropy (o/f > a-fn') f° r the most 
elongated particles in contrast to the tangential force 
anisotropy which is below the normal force anisotropy 
(aft < af n ). As a result of the increasing activation of 
tangential forces, the fraction of sliding contacts (i.e. con- 
tacts where \ft\ = (J-\f n \) grows with 77 as shown in Fig. 
l20l The contributions of side-to-side and cap-to-side con- 
tacts to force anisotropy and friction mobilization, which 
are major effects of particle shape, will be analyzed in 
section I VIII 



VI. FORCE DISTRIBUTIONS 

The force chains and spatially inhomogeneous stress 
distributions are well-known features of granular me- 
dia. A well-known observation is that a large number 
of contacts transmit very weak forces, a signature of 
the arching effect, whereas a smaller fraction of con- 
tacts carry strong force chains [3lj . Force transmis- 
sion has been investigated by experiments and numer- 
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FIG. 19: Map of radial forces for r] = 0.2 (up) and rj = 0.7 
(down) . Line thickness is proportional to the radial force. We 
represent the strong network in black and the weak network 
in red lines (see section IVlj) . The floating particles excluded 
from the force network are in white. 
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FIG. 20: Proportion of sliding contacts as a function of r\ 
averaged in the residual state. The error bars represent the 
standard deviation in the residual state. 



ical simulations for disks, ellipses and polygonal particles 
in 2D as for s phe rical, cylindrical and polyhedral parti- 
cles in 3D [ll|, [H, |H, W^h- In close correlation with 
shear strength and solid fraction, the stress transmission 
is strongly influenced by particle shape. In particular, 
one expects that elongated particle shapes will influence 
mainly the distribution of weak forces by enhancing the 
arching effect. 

The probability density function (pdf ) of radial forces 
normalized by the mean radial force (/„') is shown in 
Fig. [5l] in log-linear and log-log scales at large strains 



FIG. 21: Probability distribution function of radial forces f n i 
normalized by the average radial force (/„/ ) in log-linear (up) 
and log-log (down) plots for different values of rj. 



(the data are cumulated from several snapshots in the 
residual state) for all values of rj. As usually observed, in 
all packings the number of forces above the mean (/„/) 
falls off exponentially whereas the number of forces below 
the mean vary as a power-law: 



P(fn 



e-v(-i)H./(/.l) , /„, > (/„,), 

K H,Av) (28) 



(&yn 



fn' < (/«')> 



where a n i{rf) and (3 n '(jl) are the exponents which de- 
crease with r\ from a n (0) ~ 1.69 to a™' (0.7) ~ 0.88, and 
from (3 n '{0) ~ 0.13 to /3™'(0.7) ~ -0.49. Figure [22] shows 
the pdf P(ft') of orthoradial forces normalized by by the 
mean orthoradial force (/c) in each packing. These dis- 
tributions are also characterized by an exponential falloff 
for the forces above the average force (f t >) and a power 
law for the forces below (ft')' 



P(ft>) oc 



e-«v<"X 1 -IA'l/<IA'l>) 5 |/ t ,| > <|/ t ,|), 



, \ft'\ < 




(29) 



with the corresponding exponents ay (r/) and Pt'(v)i 
which decrease from a* (0) ~ 1.09 to a* (0.7) ~ 0.73, 
and from /3*'(Q) ~ -0.37 to /3*'(0.7) ~ -0.73. 

These distributions show clearly the larger inhomo- 
geneity of stress transmission in a granular packing com- 
posed of elongated particles. We find that (the results 
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FIG. 22: Probability distribution function of orthoradial 
forces f t i normalized by the average orthoradial force (\ft'\) 
in log-log for all values of rj. 



not shown here), as with circular particles, the contacts 
can be classified into strong and weak networks. Eval- 
uating q/p separately for each network, it is found that 
the shear stress is almost totally sustained by the strong 
contact network. 

Maps of strong and weak radial networks for radial 
forces are displayed in Fig. [19] for rj — 0.2 and rj = 0.7. 
The fraction of floating particles (less than two contacts) 
decreases linearly with r\ from 17% for r\ = to 10% 
for rj = 0.7. Hence, more particles are involved in the 
contact network for more elongated particles, but it is 
remarkable that the proportion of weak forces grows from 
60% for r\ = to 70% for r\ = 0.7. In other words, 
although the number of strong contacts decreases with r/, 
stronger force chains occur with more elongated particles. 
Although we focused here on the networks of radial force 
components, we basically obtain the same conclusions for 
the normal force components. 

We also remark that the packings are increasingly more 
inhomogenous in the sense that as particle elongation in- 
creases, the packing involves less strong force chains in 
number but with stronger forces. This decreasing force 
homogeneity in spite of increasing connectivity (fig. II 1 [) , 
means that force distributions are controlled by more 
subtle details of the microstructure than the density of 
contacts or solid fraction. As we shall see below, this is 
related to the role of various contact types in the contact 
network. 
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FIG. 23: Proportions of side-to-side (ss), cap-to-side (cs) and 
cap-to-cap (cc) contacts as a function of rj in the residual 
state. The error bars represent the standard deviation in the 
residual state. 



VII. EFFECT OF CONTACT TYPES 

Remembering that RCR particles are clumps of two 
disks with one rectangle, in this section we revisit the 
results of the previous sections in the light of the organi- 
zation of cap-to-cap (cc), cap-to-side (cs) and side-to-side 
(ss) contacts in each packing. The side-to-side or side- 
to-cap contacts do not transmit torques, but they are 
able to accommodate force lines that are usually unsus- 
tainable by cap-to-cap contacts. For this reason, it is 
worth while trying to isolate their respective roles with 
respect to the shear strength. The proportions of these 
contact types and their contributions to the structural 
anisotropy and force transmission are key quantities for 
understanding the effect of particle shape on the shear 
strength properties of granular media [llj, |35j . 

In the residual state, the proportions of different con- 
tact types are nearly constant. Fig. [231 shows the propor- 
tions fc cc , k cs and k ss of cc, cs and ss contacts averaged 
over the residual state as a function of rj. We see that 
k cc declines with 77 from 1 (for disks) to 0.2 for r\ = 0.7. 
At the same time, k cs and k ss increase from to 0.6 and 
to 0.2, respectively. Interestingly, k cs ~ k cc for 77 ~ 0.4, 
and k ss ~ k cs for 77 = 0.7. In this way, as the particle 
elongation increases, the packing passes from a contact 
network dominated by cc contacts to a contact network 
dominated by the complex contacts cs and ss. 

To identify the impact of each contact type on the 
shear strength, we proceed by additive decomposition of 
the stress tensor by considering the expression (|25p of the 
stress tensor and grouping the contacts according their 
types: 

= & cc + & cs + CTss, (30) 

where <r cc , <x cs and <r ss are obtained from the expression 
of the stress tensor Eq. ([25)) by restricting the summation 
to cc, cs and ss contacts, respectively. The correspond- 
ing stress deviators q cc , q cs and q ss are then calculated 
and normalized by the mean pressure p. Fig. [2U shows 
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FIG. 24: Shear strength (q/p)* for cs, ss and cc contacts as a 
function of rj, together with the harmonic approximation fits 
in (n,t) frame ( ) and (n',t') frame (...). 

Qcc/p, Qcs/p and q ss /p averaged in the residual state as a 
function of rj. We see clearly that q cc /p follows a trend 
opposite to that of q cs /p. For rj < 0.3, {q/p)* is dom- 
inated by cc contacts. For r\ ~ 0.3, cc and cs contacts 
participate equally to the shear stress, and for r\ > 0.3, 
the cs contacts dominate (q/p)*. Remarkably, q ss /p — 
for rj < 0.4. As we shall see below, the ss contacts par- 
ticipate to the strong force chains only in the case of the 
most elongated particles. In this way, the growth of the 
number of cs and ss contacts shown in Fig. [531 is clearly 
at the origin of a gradual consolidation of the packings 
as rj increases. 

In order to get further insight into the organization 
of different contact types, it is useful to consider partial 
connectivities P cs (c), P cc (c) and P ss (c) defined as the 
fraction of particles with exactly c contacts of cs type, 
cc type and ss type. These functions are displayed in 
Fig. US] for all our packings. Note that, by definition 
P cc {c) = P(c) for Tj — 0. We see that P cc gets narrower 
as j] increases whereas P cs and P ss get broader. It should 
be noted that, even for the most elongated particles, a 
particle can have at most two ss. This means that, the 
elongated particles tend mainly to pile up like bricks. On 
the other hand, the peak of P cs slides to the larger values 
as rj increases. For rj — 0.7, most particles have three or 
four cs contacts (for nearly 40%). 

We now consider the anisotropy of the branch vectors 
and contact forces supported by the three contact types 
at the contact and branch frames. Following the same 
procedure as for the stress tensor (see equation (|3T)|) . we 
perform an additive decomposition of the fabric and force 
tensors: 
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(31) 



where the indices refer to the partial contributions of cc, 
cs and ss contacts. The corresponding anisotropies of 
each tensor can be extracted. In principle, the principal 
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FIG. 25: Partial connectivity diagrams for all packings in the 
residual state. 

directions of these partial tensors do not coincide with 
those of the overall tensors at all stages of shearing. But, 
in practice, in the residual state, the principal directions 
coincide so that the global anisotropy of each tensor is 
the sum of its partial anisotropies: 

g 7 J ^( a c7 + a in-y + a;t 7 + a/n 7 + a/t 7 ) in (n, t) 

P [ h( a 'cy + 0/n'7 + «/n'7 + a/t'-y) in («',*')- 

(32) 

where 7 stands alternatively for {cc, cs, ss}. This decom- 
position is nicely verified by our numerical date as shown 
in FigUl 

Since the contact orientation anisotropy expressed in 
(to, t) frame and the force anisotropy expressed in (to', t') 
frame provide respectively fine descriptions of the geo- 
metrical and force organizations (see section [Vj , we re- 
strict here our analysis to the contribution of various con- 
tact types to the contact orientation anisotropy a c and 
the radial force anisotropies a t n i and a ff ■ Figure 
shows the variation of the partial contact anisotropies 
Q-ccc, o-ccs and a css due to cc, cs and ss contacts in the 
residual state as the function of rj. The anisotropy a css 
of ss contacts increases slowly with rj from to 0.18. At 
the same time, a ccc decreases and at rj = 0.7 we have 
a css = a ccc . Hence, although the ss contacts represent at 
rj = 0.7 nearly 20% of contacts, their contribution to the 
contact anisotropy remains modest and of the same order 
as cc contacts. The variation of the contact anisotropy 
a c is thus largely governed by that of a ccs . 

Figure [27] shows the partial radial force anisotropies 
a/n'cc, o-fn'cs and af n > ss , as well as the partial orthoradial 
force anisotropies af t > cc , a/t'es and a/f ss in the residual 
state as the function of rj. As for contact anisotropies, the 
cs contacts carry most of the radial and orthoradial force 
anisotropies. The ss contacts contribute modestly to the 
global force anisotropies only for 77 > 0.4. The anisotropy 
o/n'cc declines with rj, mainly due to their low number, 
and dft'cc stays nearly constant. 

A map of contact forces projected along the branch 
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FIG. 26: Partial contact orientation anisotropics a ccc , a ccs 
and a cas of cc, cs and ss contacts as the function of r\ in the 
residual state. The error bars represent the standard devia- 
tion in the residual state. 
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residual state. The error bars represent the standard devia- 
tion in the residual state. 



vectors is displayed in Fig. [55] in different colors accord- 
ing to the type of contact. For r\ — 0.7, we see that the 
network of very strong zigzag force chains is composed 
mostly of cs and ss-contacts and occasionally mediated 
by cc contacts. In contrast, for rj = 0.2, the cc contacts 
appear clearly to be correlated in the form of long chains 
across the packing rarely mediated by cs-contacts. In all 
cases, the strong force chains are mostly parallel to the 
direction of compression. 

In order to recognize quantitatively the roles of cc, cs 
and ss contacts with respect to the force network, we 
plot in Fig. |2"51 their respective proportions k cc , k cs and 
k ss alternatively for the strong and weak networks in the 
residual state ciS clS el function 77. Notice that the data 
are normalized for each network. The proportion of cc 
contacts declines rapidly in both networks as rj increases 




FIG. 28: (Color online) Snapshot of radiale forces for r\ = 0.2 
(up) and rj = 0.7 (down). Line thickness is proportional to 
the radial force. The cap-to-cap, cap-to-side and side-to-side 
contacts are in black, in red (dark gray) and in green (light 
gray). 




FIG. 29: Proportions of cap-to-cap (kcc), cap-to-side (fees) 
and side-to-side (fc B3 ) contacts in the strong (plain line) and 
week (dashed line) networks as a function of 77 in the residual 
state. The error bars represent the standard deviation in the 
residual state. 



whereas that of cs and ss contacts grow. We also remark 
that the cs contacts are slightly more numerous in the 
weak network than in the strong network. The propor- 
tions have nearly the same value in the two networks for 
cc and ss contacts. 
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VIII. CONCLUSIONS 

In this paper, we investigated the effect of particle elon- 
gation on the quasistatic behavior of sheared granular 
materials by means of Contact Dynamics simulations. 
The particle shapes are rounded-cap rectangles (RCR) 
characterized by their elongation rj defined as deviation 
from a reference circular shape, or alternatively by their 
aspect ratio. As the elongation increases from to 1, the 
particle shape varies continuously from a disk to an in- 
creasingly thin rectangle with rounded caps. The macro- 
scopic and microstructural properties of several packings 
of 13000 particles, subjected to biaxial compression, were 
analyzed as a function of r\. 

An interesting finding is that the shear strength is an 
increasing linear function of elongation, suggesting that 
the parameter ij is a "good" shape parameter for our 2D 
granular packings. In order to understand the micro- 
scopic origins of this behavior, we performed an additive 
decomposition of the stress tensor based on a harmonic 
approximation of the angular variation of average local 
branch vectors, contact normals and forces. This approx- 
imation of the shear strength in terms of texture and 
force anisotropics turns out to be in excellent agreement 
with our numerical data in the investigated range of the 
elongation parameter (j) s [0,0.7]). Given the evolution 
of various anisotropies with particle elongation, we find 
that both force and texture anisotropies contribute to the 
increase of shear strength, but the increasing mobiliza- 
tion of friction force and the associated anisotropy seem 
to be the key effect of particle elongation. In particu- 
lar the proportion of sliding contacts increases strongly 
as the particles become more elongated. This effect is 
correlated with a local nematic ordering of the particles 
which tend to be oriented perpendicular to major prin- 
cipal stress direction. This ordering is enhanced beyond 
i] = 0.4 but remains essentially of local nature. In this 
respect, the fraction of side-to-side contacts increases at 
large particle elongations. However, the force transmis- 
sion is found to be mainly guided by cap-to-side contacts, 



which represent the largest fraction of contacts for the 
most elongated particles and carry a large part of the 
shear strength. 

In contrast to shear strength, the solid fraction is not 
a monotonous function of particle elongation; It first in- 
creases with particle elongation, then declines as the par- 
ticles become more elongated. In other words, small devi- 
ation from circular shape favors the space- filling aptitude 
of the particles but beyond a characteristic elongation the 
excluded-volume effects prevail and lead to increasingly 
larger pores which cannot be filled by the particles. It is 
remarkable that the coordination number does not follow 
the solid fraction but increases with particle elongation, 
so that the packings of the most elongated particles are 
loose but well connected. 

Some features discussed in this paper can legitimately 
be attributed to the two-dimensional geometry of the 
particles. For example, rounded-cap-cylinders (sphero- 
cylinders), as three-dimensional analogs of RCR particles 
in 2D, do not undergo spontaneously a nematic order- 
ing. However, we expect that 3D packings of rounded- 
cap-cylinders behave in many ways as our 2D packings 
with increasing particle elongation. In particular, the ex- 
cluded volume effect is reinforced by particle elongation 
and it leads to a similar nonmonotonous dependance on 
the elongation as in 2D [2(J [HJ . In any case, it would 
be highly instructive to investigate 3D packings of elon- 
gated particles, along the same lines as in this paper. It is 
also obvious that more work required to assess the proper 
role of friction in 2D for RCR particles since friction mo- 
bilization seems to underlie to a large extent the shear 
strength, we are performing presently more simulations 
with lower values of friction coefficient. In the same way, 
we consider the effect of cohesion between particles with 
respect to the shear strength of packings of elongated 
particles. 

We specially thank I. Zuriguel and F. Dubois 
for fruitful discussions. This work was done as 
part of PPF CEGEO research project (www.granulo- 
science.org/CEGEO). 



[1] V. Richefeu, M. elYoussoufi, and F. Radjai, Phys. Rev. 
E 73, 051304 (2006). 

[2] V. Richefeu, E. Azema, F. Radjai, and S. Youssoufi, Pow- 
der Technology 190, 258263 (2009). 

[3] C. Nouguier-Lehon, B. Cambou, and E. Vincens, Int. J. 
Numer. Anal. Meth. Geomech 27, 1207 (2003). 

[4] S. Antony and M. Kuhn, International Journal of Solids 
and Structures 41, 5863 (2004). 

[5] N. Kruyt and L. Rothenburg, Mechanics of Materials 36, 
1157 (2004). 

[6] C. Nouguier-Lebon, E. Vincens, and B. Cambou, Inter- 
national journal of Solids and Structures 42, 6356 (2005). 

[7] F. Alonso-Marroquin, S. Luding, H. J. Herrmann, and 
I. Vardoulakis, Phys. Rev. E 71, 051304 (2005). 

[8] A. Pena, A. Lizcano, F. Alonso-Marroquin, and H. Her- 
rman, Int. J. For Numer. Anal. Meth. Geomech. 00, 1 



(2006). 

[9] C. Voivret, F. Radjai, J.-Y. Delenne, and M. S. E. Yous- 
soufi, Phys. Rev. Lett. 102, 178001 (2009). 

[10] E. Azema, F. Radjai, R. Peyroux, F. Dubois, and 
G. Saussine, Phys. Rev. E 74, 031302 (2006). 

[11] E. Azema, F. Radjai, and G. Saussine, Mechanics of Ma- 
terials 41, 721 (2009). 

[12] J. M. E. Markland, Geotechnique 31, 3,367 (1981). 

[13] Wu and Thompson, J Acoust Soc Am 108, 1046 (2000). 

[14] W. Lim and G. MacDowel, Granular Matter 7, 19 (2005). 

[15] G. Saussine, C. Cholet, P. Gautier, F. Dubois, C. Bo- 
hatier, and J. Moreau, Comput. Methods Appl. Mech. 
Eng. 195, 2841 (2006). 

[16] S. Lobo-Guerrero and L. E. Vallejo, Granular Matter 8, 
195 (2006). 

[17] M. Lu and G. McDowel, Granular Matter 9, 69 (2007). 



17 



[18] A. Tkachenko and T. A. Witten, Phys Rev E 60, 687 
(1999). 

[19] A. Donev, F. Stillinger, P. Chaikin, and S. Torquato, 

Phys Rev Lett. 92, 255506 (2004). 
[20] A. Donev, I. Cisse, D. Sachs, E. Variano, F. Stillinger, 

R. Connelly, S. Torquato, and P. Chaikin, Science 303, 

990 (2004). 

[21] W. Man, A. Donev, F. Stillinger, M. Sullivan, W. Russel, 

D. Heeger, S.Inati, S. Torquato, and P. Chaikin, Phys 

Rev Lett pp. 198001-1, 198001-4 (2005). 
[22] A. Donev, R. Connelly, F. Stillinger, and S. Torquato, 

Phys Rev E 75, 051304 (2007). 
[23] G. Yatsenko and K. Schweizer, Phys Rev E 76, 041506 

(2007). 

[24] H. Ouadfel and L. Rothenburg, Mechanics of Materials 
33, 201 (2001). 

[25] I. Zuriguel, T. Mullin, and J. Rotter, Phys. Rev. Lett. 
98, 028001 (2007). 

[26] R. C. Hidalgo, I. Zuriguel, D. Maza, and I. Pagonabar- 
raga, Phys. Rev. Lett. 103, 118001 (2009). 

[27] A. Kyrylyuk, A. Wouterse, and A. Philipse, AIP Proc. 
1145, 211 (2009). 

[28] N. P. Kruyt and L. Rothenburg, ASME Journal of Ap- 
plied Mechanics 118, 706 (1996). 

[29] R. J. Bathurst and L. Rothenburg, J. Appl. Mech. 55, 
17 (1988). 

[30] L. Rothenburg and R. J. Bathurst, Geotechnique 39, 601 
(1989). 

[31] F. Radjai, D. E. Wolf, M. Jean, and J. Moreau, Phys. 
Rev. Letter 80, 61 (1998). 

[32] A. Mirghasemi, L. Rothenburg, and E. Maryas, Geotech- 
nique 52, N 3, 209 (2002). 

[33] H. Troadec, F. Radjai, S. Roux, and J.-C. Charmet, 
Phys. Rev. E 66, 041305 (2002). 

[34] F. Radjai, in 6th International Conference on the Mi- 
cromechanics of Granular Media, JUL 13-17, 2009 
Golden, CO POWDERS AND GRAINS 2009 (AIP Con- 
ference Proceedings, 2009), vol. 1145, pp. 35-42. 

[35] E. Azema, F. Radjai, R. Peyroux, and G. Saussine, Phys. 
Rev. E 76, 011301 (2007). 

[36] S. Sacanna, L. Rossi, A. Wouterse, and A. P. Philipse, J. 
Phys.: Condens. Matter 19 (2007) 376108 376108, 16p 
(2007). 

[37] S. Williams and A. Philipse, Phys. Rev. E 67, 051301 
(2003). 

[38] F. Fraige, P. . Langston, and G. Chen, Powder Technol- 
ogy 186, 224 (2008). 

[39] P. Langston, M. Al-Awamleh, F. Fraige, and B. Asmar, 
Chemical Engineering Science 59, 425 (2004). 

[40] M. Jean and J. J. Moreau, in Proceedings of Contact 
Mechanics International Symposium (Presses Polytech- 
niques et Universitaires Romandes, Lausanne, Switzer- 
land, 1992), pp. 31-48. 

[41] J. Moreau, Eur. J. Mech. A/Solids 13, 93 (1994). 

[42] M. Jean, Computer Methods in Applied Mechanic and 
Engineering 177, 235 (1999). 

[43] J. Moreau, in Novel approaches in civil engineering, 
edited by M. Fremond and F. Maceri (Springer- Verlag, 
2004), no. 14 in Lecture Notes in Applied and Computa- 
tional Mechanics, pp. 1-46. 

[44] F. Radjai and S. Roux, in He Congres Francais de 
Mecanique. Toulouse (1999). 

[45] F. Dubois and M. Jean, in Actes du sixieme col- 
loque national en calcul des structures - CSMA- 



AFM-L MS - (2003), vol. 1, pp. 111-118, URL 
https : / / subver . lmgc .univ-montp2 . f r/trac_LMGC90v2/ 

[46] M. Renouf and P. Alart, Comput. Methods Appl. Mech. 
Engrg 194, 2019 (2005). 

[47] F. Radjai, Physics of dry granular media (Kluwer Aca- 
demic Publishers (Dordrecht/Boston/London), 1997), 
chap. Multicontacts dynamics, p. 305. 

[48] F. Radjai and E. Azema, Eur. J. Env. Civil Engineering 
13, 204 (2009). 

[49] F. Radjai and V. Richefeu, Mechanics of Materials 41, 
715 (2009). 

[50] R. Folk, Petrology of Sedimentary Rocks (Hemphill Pub- 
lishing Company, Austin Texas 78703, 1974). 

[51] C. Cholet, G. Saussine, P. Gautier, F. Dubois, C. Bo- 
hatier, G. Combe, and K. Sab, in World Congress on 
Railway Research (WCRR) (2003). 

[52] E. Azema, F. Radjai, R. Peyroux, V. Richefeu, and 
G. Saussine, Eur. Phys. J. E 26, 327 (2008). 

[53] I. Bratberg, F. Radjai, and A. Hansen, Phys. Rev. E 66, 
031303 (2002). 

[54] A. Taboada, K. J. Chang, F. Radjai, and F. Bouchette, 
Journal Of Geophysical Research 110, 1 (2005). 

[55] C. Voivret, F. Radjai, J.-Y. Delenne, and M. S. E. Yous- 
soufi, Phys Rev E 76, 021301 (2007). 

[56] C. Voivret, Ph.D. thesis, Univer- 
sity Montpellier II (2008), URL 

|http : //tel . archives- ouvertes . f r/tel -00372125 _vl/ 1 

[57] GDR-MiDi, Eur. Phys. J. E 14, 341 (2004). 

[58] J. J. Moreau, in Friction, Arching, Contact Dynamics, 
edited by D. E. Wolf and P. Grassberger (World Scien- 
tific, Singapore, 1997), pp. 233-247. 

[59] L. Staron and F. Radjai, Phys. Rev. E 72, 1 (2005). 

[60] J. Mitchell and K. Soga, Fundamentals of Soil Behavior 
(Wiley, New- York, NY, 2005). 

[61] J. Blouwolff and S. Fraden, Europhys. Lett. 76, 1095 
(2006). 

[62] A. Wouterse, S. Williams, and A. Philipse, J. Phys.: Con- 
dens. Matter 19 406215, 14 (2007). 

[63] D. Wood, Soil behaviour and critical state soil mechan- 
ics (Cambridge University Press, Cambridge, England, 
1990). 

[64] F. Radjai and S. Roux, in The Physics of Granular Me- 
dia, edited by H. Hinrichsen and D. E. Wolf (Wiley- VCH, 
Weinheim, 2004), pp. 165-186. 

[65] D. Taylor, Fundamentals of soils mechanics (Wiley, New 
York,' 1948). 

[66] P. Vermeer, in Constitutive relation for soils (Rotterdam, 
A. A Balkema, 1984), pp. 175-197. 

[67] P. Vermeer and R. Borst, Heron 29, 1 (1984). 

[68] A. Taboada, N. Estrada, and F. Radjai, Phys. Rev. Lett. 
97, 098302 (2006). 

[69] L. Rothenburg, R. Bathurst, and A. Berlin, in Powder 
and Grains, edited by Thornton (Balkema, Rotterdam, 
1993), pp. 147-153. 

[70] K. A. Reddy, V. Kumaran, and J. Talbot, Phys Rev E 
80, 031304 (2009). 

[71] I. Zuriguel and T. Mullin, Proc. R. Soc. A 464, 99 (2008). 

[72] J.-N. Roux, Phys Rev E 61, 6802 (2000). 

[73] C. F. Moukarzel, Phys Rev Letter 81, 1634 (1998). 

[74] S. F. Edwards, Physica A 249, 226 (1998). 

[75] S. Alexander, Phys. Rep. 296, 65 (1998). 

[76] R. Guises, J. Xiang, J. -P. Latham, and A. Munjiza, Gran- 
ular Matter 11, 281 (2009). 



18 



[77] M. Mailman, C. F. Schreck, C. OOHern, and [85 

B. Chakraborty, Phys. Rev. Letters 102, 255501 (2009). 
[78] M. Satake, in Proceedings of the IUTAM symposium 

on deformation and failure of granular materials, Delft, 
edited by P. A. Vermeer and H. J. Luger (A. A. Balkema, 
Amsterdam, 1982), pp. 63-68. [87 

[79] M. Oda and K. Iwashita, eds., Mechanics of Granular 

Materials (A. A. Balkema, Rotterdam, 1999). [ 

[80] B. Cambou, P. Dubujet, and C. Nouguier-Lehon, Me- 
chanics of Materials 36, 1185 (2004). [89 

[81] A. Pena, R. Garcia-Rojo, and H. Herrmann, Granular 

Matter 9, 279 (2007). ' [90 

[82] F. Radjai, H. Troadec, and S. Roux, in Granular Materi- 
als: Fundamentals and Applications, edited by S. Antony, [91 
W. Hoyle, and Y. Ding (RS.C, Cambridge, 2004), pp. [92 
157-184. 

[83] M. Oda, J. Koshini, and S. Nemat-Nasser, Geotechnique [93 
30, 479 (1980). 

[84] B. Cambou, in Powders and Grains 93, edited by [94 

C. Thornton (A. A. Balkema, Amsterdam, 1993), pp. 
73-86. 



L. Landau and L. Lifshit, Statistical physics : course of 
theoretical physics (Oxford Pergamon Press 1-2, 1959). 

C. Liu, S. R. Nagel, D. A. Schecter, S. N. Coppersmith, 
S. Majumdar, O. Narayan, and T. A. Witten, Science 
269, 513 (1995). 

D. M. Mueth, H. M. Jaeger, and S. R. Nagel, Phys. Rev. 

E. 57, 3164 (1998). 

F. Radjai, M. Jean, J. Moreau, and S. Roux, Phys. Rev. 
Letter 77, 274 (1996). 

G. Lovol, K. Maloy, and E. Flekkoy, Phys. Rev. E 60, 
5872 (1999). 

S. G. Bardenhagen, J. U. Brackbill, and D. Sulsky, Phys. 
Rev. E 62, 3882 (2000). 

S. J. Antony, Phys Rev E 63, 011302 (2001). 

L. E. Silbert, G. S. Grest, and J. W. Landry, Phys. Rev. 

E 66, 1 (2002). 

T. S. Majmudar and R. P. Behringer, Nature 435, 1079 
(2005). 

P. T. Metzger, Phys Rev E 77, 011307 (2008). 



